// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

#include <iostream>
#include <ostream>
#include <sstream>
#include <vector>
#include <algorithm>
#include <complex>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <algorithm>
#include <iomanip>
#include <vector>
#include <list>
#include <set>
#include <cstdarg>
#include <stdexcept>

#include <Eigen/Eigen>

#include <cmdioutil.h>

void load_hamiltonian_old(std::string fn,std::vector<std::pair<int,int> > & range,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & h,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & s,double &ef_m)
{
    std::complex<double> ii(0,1.0);
    
    std::vector<std::string> lines;
    io::load(fn,"#",lines);

    ef_m = 0.0;
    std::vector<std::string> efermi_block;
    io::parse_block("<fermi-energy>","</fermi-energy>",lines,efermi_block);
    for(std::vector<std::string>::iterator it = efermi_block.begin(); it != efermi_block.end(); it++)
    {
        std::istringstream iss(*it);
        if (!(iss >> ef_m)) ef_m = 0.0;
    }

    
    std::vector<std::string> orbital_block;
    io::parse_block("<index-orbital>","</index-orbital>",lines,orbital_block);
    
    range.clear();
    int index,orbital;
    int orbital_index = 0;
    for(std::vector<std::string>::iterator it = orbital_block.begin(); it != orbital_block.end(); it++)
    if(util::match("% %",*it,index,orbital))
    {
            std::pair<int,int> p(orbital_index,orbital_index + orbital);
            orbital_index += orbital;
            range.push_back(p);
    }
    
    std::vector<std::string> hamiltonian_block;
    io::parse_block("<hamiltonian-overlap>","</hamiltonian-overlap>",lines,hamiltonian_block);
    
    int n = orbital_index;
    h.resize(n,n); h.setZero(); s = h;
    double res,ims,reh,imh;
    int i,j;
    for(std::vector<std::string>::iterator it = hamiltonian_block.begin(); it != hamiltonian_block.end(); it++)
    if(util::match("% % % % % %",*it,i,j,res,ims,reh,imh))
    {
        h(i-1,j-1)=reh+ii*imh;
        s(i-1,j-1)=res+ii*ims;
    }

}


void load_hamiltonian(std::string fn,std::vector<std::pair<int,int> > & range,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & h,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & s,double &ef_m)
{
    std::complex<double> ii(0,1.0);
    
    std::vector<std::string> lines;
    io::load(fn,"#",lines);

    ef_m = 0.0;
    std::vector<std::string> efermi_block;
    io::parse_block("<fermi-energy>","</fermi-energy>",lines,efermi_block);
    for(std::vector<std::string>::iterator it = efermi_block.begin(); it != efermi_block.end(); it++)
    {
        std::istringstream iss(*it);
        if (!(iss >> ef_m)) throw std::runtime_error("Not able to parse fermi energy block in the input file!");
        break;
    }

    std::vector<std::string> orbital_block;
    if(!io::parse_block("<index-orbital>","</index-orbital>",lines,orbital_block)) throw std::runtime_error("No or inproper index-orbital block in the input file!");
    if(!orbital_block.size()) throw std::runtime_error("No or inproper index-orbital block in the input file!");

    
    range.clear();
    int index,orbital;
    int orbital_index = 0;
    for(std::vector<std::string>::iterator it = orbital_block.begin(); it != orbital_block.end(); it++)
    if(util::match("% %",*it,index,orbital))
    {
            std::pair<int,int> p(orbital_index,orbital_index + orbital);
            orbital_index += orbital;
            range.push_back(p);
    }
    int n = orbital_index;
    h.resize(n,n); h.setZero(); s = h; s.setIdentity();
    double re,im;
    int i,j;

    std::vector<std::string> overlap_block;
    if(!io::parse_block("<overlap>","</overlap>",lines,overlap_block)) throw std::runtime_error("There is no overlap matrix block or it is inproper in the input file!");
    for(std::vector<std::string>::iterator it = overlap_block.begin(); it != overlap_block.end(); it++)
    if(util::match("% % % % ",*it,i,j,re,im)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) s(i-1,j-1)=re+ii*im; else throw std::runtime_error("Inconsistent overlap matrix index with index orbital block!");  }
//    else if(util::match("% % % ",*it,i,j,re)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) s(i-1,j-1)=re+ii*0.0; else throw std::runtime_error("Inconsistent overlap matrix index with index orbital block!"); }
    
    std::vector<std::string> hamiltonian_block;
    if(!io::parse_block("<hamiltonian>","</hamiltonian>",lines,hamiltonian_block)) throw std::runtime_error("There is no hamiltonian block or it is inproper in the input file!");
    for(std::vector<std::string>::iterator it = hamiltonian_block.begin(); it != hamiltonian_block.end(); it++)
    if(util::match("% % % % ",*it,i,j,re,im)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) h(i-1,j-1)=re+ii*im; else throw std::runtime_error("Inconsistent hamiltonian index with index orbital block!"); }
//    else if(util::match("% % % ",*it,i,j,re)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) h(i-1,j-1)=re+ii*0.0; else throw std::runtime_error("Inconsistent hamiltonian index with index orbital block!"); }
    if(!hamiltonian_block.size()) throw std::runtime_error("There is no hamiltonian block or it is inproper in the input file!");

}


void load_hamiltonian_gollum(std::string fn,std::vector<std::pair<int,int> > & range,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & h,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & s,double &ef_m)
{
    std::ifstream is(fn.c_str());
    std::complex<double> ii(0,1.0);
    std::string dump;
    is >> dump >> dump >> dump >> dump >> dump >> dump >> dump;
    is >> dump >> dump >> dump >> dump >> dump >> dump >> ef_m;
    is >> dump >> dump >> dump >> dump >> dump >> dump;
    int rows,cols; 
    is >> dump >> dump >> rows;
    is >> dump >> dump >> cols;
    int index;
    int pindex = 1;
    int pi = 0;
    int n = rows;    
    for(int j = 0; j < n; j++)
    {
        is >> index >> dump >> dump >> dump >> dump;
        if(index != pindex)
        {
            std::pair<int,int> p(pi,j);
            range.push_back(p);
            pindex = index;
            pi = j;
        }
    }
    std::pair<int,int> p(pi,n);
    range.push_back(p);

    is >> dump >> dump >> dump >> dump >> dump >> dump;
    //int rows,cols;
    is >> dump >> dump >> rows;
    is >> dump >> dump >> cols;
    for(int j = 0; j < rows; j++)
        is >> dump >> dump >> dump;

    is >> dump >> dump >> dump >> dump >> dump >> dump;
    //int rows,cols;
    is >> dump >> dump >> rows;
    is >> dump >> dump >> cols;
    h.resize(n,n); h.setZero(); s = h;
    for(int k = 0; k < rows; k++)
    {
        int i,j;
        double res,ims,reh,imh;
        is >> dump >> i >> j >> res >> ims >> reh >> imh;
        h(i-1,j-1)=reh+ii*imh;
        s(i-1,j-1)=res+ii*ims;

    }
    is.close();
    
}


#ifdef DEBUG
void out_complex_matrix(Eigen::Matrix< std::complex<double> , Eigen::Dynamic, Eigen::Dynamic >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        //os.flags(std::ios::scientific);
        os.flags(std::ios::fixed);
        if(comment)
        {
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "Re Im ... Re Im" << ' '; 
                os << std::endl;
        }
        for(int i = 0; i < m.rows(); i++)
        {
	        for(int j = 0; j < m.cols(); j++)
                {
                        os << std::right << std::setw(w+4) << std::setprecision(p) << real(m(i,j)) << ' '; 
                        os << std::right << std::setw(w) << std::setprecision(p) << imag(m(i,j)) << ' '; 
                }
                os << std::endl;                        
        }                
        os << std::endl;                        
}
#endif

#ifdef DEBUG
void out_complex_vector(Eigen::Matrix< std::complex<double> , Eigen::Dynamic, 1 >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        os.flags(std::ios::scientific);
        if(comment)
        {
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "Re Im" << ' '; 
                os << std::endl;
        }
        for(int i = 0; i < m.rows(); i++)
        {
                        os << std::right << std::setw(w+4) << std::setprecision(p) << real(m(i)) << ' '; 
                        os << std::right << std::setw(w) << std::setprecision(p) << imag(m(i)) << ' '; 
                os << std::endl;                        
        }                
        os << std::endl;                        
}
#endif

#ifdef DEBUG
void out_real_matrix(Eigen::Matrix< double , Eigen::Dynamic, Eigen::Dynamic >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        os.flags(std::ios::scientific);
        for(int i = 0; i < m.rows(); i++)
        {
	        for(int j = 0; j < m.cols(); j++)
                {
                        os << std::right << std::setw(w+4) << std::setprecision(p) << m(i,j) << ' '; 
                }
                os << std::endl;                        
        }                
        os << std::endl;                        
}
#endif


#ifdef DEBUG
void out_ssd(double e, Eigen::Matrix< std::complex<double> , Eigen::Dynamic, Eigen::Dynamic > & s,std::ostream & os,bool comment=true, int w = 12,int p = 5)
{
        os.flags(std::ios::scientific);
        if(comment)
        {
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "E-Ef [eV]" << ' '; 
                os << std::left << std::setw(w) << std::setprecision(p) << "maximum norm of S*adj(S)-I " << ' '; 
                os << std::endl;
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "Re Im ... Re Im" << ' '; 
                os << std::endl;
        }
                
        os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << e << ' ';
        if( (s.rows() != 0) && (s.cols() != 0) )
        {
            Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > m;
            m = s*s.adjoint();

            os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << ((m-Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic >::Identity(m.rows(),m.cols())).array().abs()).maxCoeff() << ' ';
            os << std::endl;                
            for(int i = 0; i < m.rows(); i++)
            {
    	            for(int j = 0; j < m.cols(); j++)
                    {
                            os << std::right << std::setw(w+4) << std::setprecision(p) << std::fixed << real(m(i,j)) << ' '; 
                            os << std::right << std::setw(w) << std::setprecision(p) << std::fixed << imag(m(i,j)) << ' '; 
                    }
                    os << std::endl;                        
            }                
        }
        else os << std::endl;                        
        os << std::endl;                        
}
#endif

template<typename scattering>
void out_coef(double e, scattering & s,std::ostream & os,bool comment = false,bool advanced = false,int w = 18,int p = 10)
{
        //os.flags(std::ios::scientific);
        os.flags(std::ios::fixed);
        if(comment)
        {       
                std::ostringstream oss;
                oss << "E-Ef(=" << std::setprecision(2)<< std::fixed << e <<")";
                os << '#';
                os << std::right << std::setw(w-1) << std::setprecision(p) << oss.str() << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Tran a->b" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Tran b->a" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Ref a->a" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Ref b->b" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "a oc" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "b oc" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Max tt' a->b" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Max tt' b->a" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "max abs S*S'-I" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "eiva tt' a->b" << ' '; 
                os << std::endl;
        }
        else
        {   
            os << std::right << std::setw(w) << std::setprecision(p) << e << ' ';
            if(s.ano*s.bno !=0)
            {
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefba << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefab << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefaa << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefbb << ' ';
            } 
            else
            {
                os << std::right << std::setw(w) << std::setprecision(p) << 0.0 << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << 0.0 << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << (double)s.ano << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << (double)s.bno << ' ';
                s.checknorm = 0.0;
            }
            os << std::right << std::setw(w) << std::setprecision(p) << s.ano << ' '; 
            os << std::right << std::setw(w) << std::setprecision(p) << s.bno << ' '; 
            if(s.ano*s.bno !=0)
            {
                os << std::right << std::setw(w) << std::setprecision(p) << s.maxcoefba << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.maxcoefab << ' '; 
            } 
            else
            {
                os << std::right << std::setw(w) << std::setprecision(p) << 0.0 << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << 0.0 << ' '; 
                s.checknorm = 0.0;
            }
/*
            os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << s.checknorm << ' '; 
            for(int i = s.ttabe.rows()-1; i >= 0; i--)
            os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << s.ttabe(i).real() << ' ';  
            os << std::endl;                
            
            if(advanced)
            {
                os << std::endl;
                os << "advanced: scatterer wv for scatterer channles (experiment 2010jul16)" << std::endl;
                for(int i = s.ttabe.rows()-1; i >= 0; i--)
                os << std::right << std::setw(33) << std::setprecision(10) << std::scientific << s.ttabe(i).real() << ' ';  
                os << std::endl;                

                for(int i = 0; i < s.yml.rows(); i++)
                {
                    for(int j = s.yml.cols()-1; j >= 0; j--)
                    {
                        os << std::right << std::setw(15) << std::setprecision(5) << std::scientific << s.yml(i,j).real() << ' ';  
                        os << std::right << std::setw(17) << std::setprecision(5) << std::scientific << s.yml(i,j).imag() << ' ';  
                    }
                    os << std::endl;
                }
                os << std::endl;
                os << std::endl;
            }
*/

        }
        
        

}


void get_orbital_range(const std::vector<std::pair<int,int> > & range,const util::pair_expression<int,int,':'> &expr, int &io,int &no)
{
        int index1,index2,n;
        n = range.size();
        index1 = std::get<0>(expr);
        index2 = std::get<1>(expr);
        if(index1*index2 == 0) throw std::runtime_error("Index cannot be zero.");
        int i1,i2;
        i1 = index1 < 0 ? index1 + n : index1 - 1;
        i2 = index2 < 0 ? index2 + n : index2 - 1;
        if(i1 > i2) std::swap(i1,i2);                    
    
        io = range[i1].first;
        no = range[i2].second-io;
}


template<typename lead,typename scatterer>
void get_amb(int argc, char* argv[],lead & a,scatterer & m,lead & b,double & emin_p, double & emax_p,double & ef_p,int & nume_p,double &tol_p,bool &advanced)
{
        cmd::arguments arg(argc,argv);
        cmd::switcharg<false> help(arg,"-h","--help", "print help message");
        cmd::switcharg<false> info(arg,"-i","--info", "print info/example message");

        cmd::vararg<double> emin(arg,"-emin","--emin","min of e-range (ef will be set to zero), <arg>=-1.0",-1.0);
        cmd::vararg<double> emax(arg,"-emax","--emax","max of e-range (ef will be set to zero), <arg>=1.0",1.0);
        //cmd::vararg<double> ewindow(arg,"-ew","--ewindow","e-range  window (ef will be set to zero), <arg>=1.0",0.0);
        cmd::vararg<double> ef(arg,"-ef","--ef","fermi energy, <arg>=0.0 or takes it from the hamiltonian file if there is",0.0);
        cmd::vararg<int> nume(arg,"-n","--nume","number of energy intervals, <arg>=1",1);
        
        
        cmd::vararg<util::pair_expression<int,int,':'> > ha(arg,"-ha","--ha","lead-a uc pair atom index range, 1:3 or -1:-4","");
        cmd::vararg<util::pair_expression<int,int,':'> > hb(arg,"-hb","--hb","lead-b uc pair atom index range, 1:4 or -1:30","");
        cmd::vararg<util::pair_expression<int,int,':'> > ah(arg,"-ah","--ah","lead-a uc pair atom index range (h taken from opposite side)","");
        cmd::vararg<util::pair_expression<int,int,':'> > bh(arg,"-bh","--bh","lead-b uc pair atom index range (h taken from opposite side)","");
        cmd::varargxx<util::pair_expression<int,int,':'>,std::string> la(arg,"-la","--lead-a","replace i:j h range and lead-a ham. file","","");
        cmd::varargxx<util::pair_expression<int,int,':'>,std::string> lb(arg,"-lb","--lead-b","replace i:j h range lead-b ham. file","","");

        cmd::vararg<double> tol(arg,"-tol","--tol","channel separation tolerance def.: 1e-8",1.0e-8);
        cmd::switcharg<false> rand(arg,"-r","--rand", "default diagonal randomizer");        
        cmd::vararg<double> randa(arg,"-ra","--randa","randomize lead-a, to avoid singular coupling, <arg>=0.0 ",0.0);
        cmd::vararg<double> randb(arg,"-rb","--randb","randomize lead-b, to avoid singular coupling, <arg>=0.0",0.0);


        cmd::varargxx<util::pair_expression<int,int,':'>,double> homo(arg,"-homo","--homo","homo level shift of the selected range","",0.0);
        cmd::varargxx<util::pair_expression<int,int,':'>,double> lumo(arg,"-lumo","--lumo","lumo level shift of the selected range","",0.0);

        
        cmd::switcharg<false> gollum(arg,"-g","--gollum", "hamiltonian file in gollum format");
        cmd::switcharg<false> symmetrize(arg,"-s","--sym", "symetrize hamiltonian file before calculation");
        cmd::switcharg<false> adv(arg,"-a","--adv", "Advanced analysis and outputs");

        cmd::posarg<std::string> ham_fn(arg,"hamiltonain input file name", "");

        if(*help)
        {
                arg.print_help(std::cout,"transportab");
                std::cout << "examples:" << std::endl;
                std::cout << "transportab hamiltonian.dat -ha 3:4 -hb -4:-3 -n 100" << std::endl;
                std::cout << "transportab hamiltonian.dat -ha 3:4 -hb -4:-3 -n 100 -r" << std::endl;
                std::cout << "version 0.9b" << std::endl;
                exit(EXIT_SUCCESS);
        }
        if(*info)
        {
                std::cout << "version 0.9b" << std::endl;
                std::cout << "Copyright (c) 2010 David Zsolt Manrique (david.zsolt.manrique@gmail.com)"  << std::endl;
                
                //return EXIT_SUCCESS;
                exit(EXIT_SUCCESS);
        
        }                                                                                                                
        
        if(!~ham_fn) throw std::runtime_error("Give hamiltonian file!");
        
        std::vector<std::pair<int,int> > range,arange,brange;
        Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> hfull,sfull;

        double ef_m =0.0;
        if(!*gollum) load_hamiltonian(*ham_fn,range,hfull,sfull,ef_m);
        else load_hamiltonian_gollum(*ham_fn,range,hfull,sfull,ef_m);
        
        if(*symmetrize) 
        { 
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> tmp;
            tmp = hfull.adjoint();
            hfull += tmp;
            tmp = sfull.adjoint();
            sfull += tmp;
            hfull *= 0.5;  sfull *= 0.5; 
        }
        else
        {
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> tmp;
            tmp = hfull.adjoint();
            tmp -= hfull;
            if(tmp.norm() > 1e-6 ) { std::ostringstream ss; ss << tmp.norm(); throw std::runtime_error("Hamiltonian is not hermitian! Antisym norm: " + ss.str()); }
            tmp = sfull.adjoint();
            tmp -= sfull;
            if(tmp.norm() > 1e-6 ) { std::ostringstream ss; ss << tmp.norm(); throw std::runtime_error("Overlap is not hermitian! Antisym norm: " + ss.str()); }
            /*check if gramm matrix e?*/

        }

        
        
        

        if(~homo) 
        { 
            int homo_i,homo_n;
            get_orbital_range(range,homo.var1,homo_i,homo_n); 
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> homos,homoh, ish,evec;
            Eigen::VectorXcd eval;
            homoh = hfull.block(homo_i,homo_i,homo_n,homo_n);
            homos = sfull.block(homo_i,homo_i,homo_n,homo_n);
            ish = homos.inverse()*homoh;
            Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(ish);  
            evec = eigen_solver.eigenvectors();                                             
            eval = eigen_solver.eigenvalues();                                                               
            
            for(int i = 0; i < eval.rows(); i++) 
            if(eval(i).real()<ef_m) eval(i)+=homo.var2; // na melyik fermi a command lineos vagy a filebol !! ezt konzisztense kell majd tenni
            
            
            hfull.block(homo_i,homo_i,homo_n,homo_n) = homos * evec * eval.asDiagonal() * evec.adjoint();
        }

        if(~lumo) 
        { 
            int lumo_i,lumo_n;
            get_orbital_range(range,lumo.var1,lumo_i,lumo_n); 
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> lumos,lumoh, ish,evec;
            Eigen::VectorXcd eval;
            lumoh = hfull.block(lumo_i,lumo_i,lumo_n,lumo_n);
            lumos = sfull.block(lumo_i,lumo_i,lumo_n,lumo_n);
            ish = lumos.inverse()*lumoh;
            Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(ish);  
            evec = eigen_solver.eigenvectors();                                             
            eval = eigen_solver.eigenvalues();                                                               
            
            for(int i = 0; i < eval.rows(); i++) 
            if(eval(i).real()>ef_m) eval(i)+=lumo.var2; // na melyik fermi a command lineos vagy a filebol !! ezt konzisztense kell majd tenni
            
            
            hfull.block(lumo_i,lumo_i,lumo_n,lumo_n) = lumos * evec * eval.asDiagonal() * evec.adjoint();
        }

        
        int nm=hfull.rows();
        int im=0;

        std::vector<std::pair<int,int> > larange,lbrange;
        Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> lah,las,lbh,lbs;
        double ef_la=0.0;
        double ef_lb=0.0;
        if(~la)
        {
            if(!*gollum) load_hamiltonian(la.var2,larange,lah,las,ef_la);
            else load_hamiltonian_gollum(la.var2,larange,lah,las,ef_la);
            double def = ef_m - ef_la;
            lah += def * las;
        }
        if(~lb)
        {
            if(!*gollum) load_hamiltonian(lb.var2,lbrange,lbh,lbs,ef_lb);
            else  load_hamiltonian_gollum(lb.var2,lbrange,lbh,lbs,ef_lb);
            double def = ef_m - ef_lb;
            lbh += def * lbs;
        }

        //////////////////////////////////////////////////////////////////
    
        int ia2,na2;
        int ia,na;
        if(~ha) { get_orbital_range(range,*ha,ia2,na2); ia=ia2; na=na2/2; }
        else if(~ah) { get_orbital_range(range,*ah,ia2,na2); na=na2/2; ia=ia2+na; }
        else if(~la) { get_orbital_range(range,la.var1,ia,na); }
        else throw std::runtime_error("Give -ha -ah range or -la range lead ham. file!");

        int ib2,nb2;
        int ib,nb;
        if(~hb) { get_orbital_range(range,*hb,ib2,nb2);  nb=nb2/2;ib=ib2+nb; }
        else if(~bh) { get_orbital_range(range,*bh,ib2,nb2); nb=nb2/2; ib=ib2; }
        else if(~lb) { get_orbital_range(range,lb.var1,ib,nb); }
        else throw std::runtime_error("Give -hb -bh range or -lb range lead ham. file!");
        
        ////////////////////////////////////////////
        
        //std::cout << ia << "  " << na << std::endl;
        //std::cout << ib << "  " << nb << std::endl;        
        
        if(~la)
        {
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> la_h,la_s,a_h,a_s;

            a_h = hfull.block(ia,ia,na,na);
            a_s = sfull.block(ia,ia,na,na);
        
            la_h = lah.block(na,na,na,na);
            la_s = las.block(na,na,na,na);
            
            double de = (la_h - a_h).diagonal().array().sum().real() / double(na);
            lah -= de * las;

            a.h = lah.block(na,na,na,na);
            a.s = las.block(na,na,na,na);
            a.v = lbh.block(0,na,na,na);
            a.t = lbs.block(0,na,na,na);
            
            im=ia+na;
        }
        else if(~ha)
        {
            a.h = hfull.block(ia,ia,na,na);
            a.s = sfull.block(ia,ia,na,na);
            a.v = hfull.block(ia,ia+na,na,na);
            a.t = sfull.block(ia,ia+na,na,na);
            im=ia+na;
        }
        else if(~ah)
        {
            a.h = hfull.block(ia,ia,na,na);
            a.s = sfull.block(ia,ia,na,na);
            a.v = hfull.block(ia-na,ia,na,na);
            a.t = sfull.block(ia-na,ia,na,na);
            im=ia+na;
        }


        if(~lb)
        {
        
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> lb_h,lb_s,b_h,b_s;

            b_h = hfull.block(ib,ib,nb,nb);
            b_s = sfull.block(ib,ib,nb,nb);

            lb_h = lbh.block(nb,nb,nb,nb);
            lb_s = lbs.block(0,nb,nb,nb);
        
            double de = (lb_h - b_h).diagonal().array().sum().real() / double(nb);
            lbh -= de * lbs;
        
            b.h = lbh.block(nb,nb,nb,nb);
            b.s = lbs.block(nb,nb,nb,nb);
            b.v = lbh.block(0,nb,nb,nb);
            b.t = lbs.block(0,nb,nb,nb);

            nm=ib-im;

        }
        else if(~hb)
        {
            b.h = hfull.block(ib,ib,nb,nb);
            b.s = sfull.block(ib,ib,nb,nb);
            b.v = hfull.block(ib-nb,ib,nb,nb);
            b.t = sfull.block(ib-nb,ib,nb,nb);
            nm=ib-im;
        }
        else if(~bh)
        {
            b.h = hfull.block(ib,ib,nb,nb);
            b.s = sfull.block(ib,ib,nb,nb);
            b.v = hfull.block(ib,ib+nb,nb,nb);
            b.t = sfull.block(ib,ib+nb,nb,nb);
            nm=ib-im;
        }

        
        m.h = hfull.block(im,im,nm,nm);
        m.s = sfull.block(im,im,nm,nm);
        
        m.va = hfull.block(im-na,im,na,nm);
        m.ta = sfull.block(im-na,im,na,nm);
        
        m.vb = hfull.block(im,im+nm,nm,nb);
        m.tb = sfull.block(im,im+nm,nm,nb);


        if(*rand)
        {

            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> r = a.v;
            r.setRandom();
            a.v += 0.0001 * r.real().cast<std::complex<double> > ();
            r = b.v;
            r.setRandom();
            b.v += 0.0001 * r.real().cast<std::complex<double> > ();
/*
            Eigen::VectorXcd r;
            r.resize(a.v.rows());
            r.setRandom();
            a.v += 0.001 * r.real().cast<std::complex<double> > ().asDiagonal();
            
            r.resize(b.v.rows());
            r.setRandom();
            b.v += 0.001 * r.real().cast<std::complex<double> > ().asDiagonal();
*/

        }
        
        if(~randa)
        {
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> r = a.v;
            r.setRandom();
            a.v += *randa * r.real().cast<std::complex<double> > ();
        }

        if(~randb)
        {
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> r = b.v;
            r.setRandom();
            b.v += *randb * r.real().cast<std::complex<double> > ();
        }

        tol_p = *tol;
        emin_p = *emin;
        emax_p = *emax;
        if(~ef) ef_p = *ef; else ef_p = ef_m;
        nume_p = *nume;
        advanced = *adv;
                
#ifdef DEBUG
	std::cout << "hfull:" << std::endl;
        out_complex_matrix(hfull,std::cout);
	std::cout << "sfull:" << std::endl;
        out_complex_matrix(sfull,std::cout);


	std::cout << "ha:" << std::endl;
        out_complex_matrix(a.h,std::cout);
	std::cout << "sa:" << std::endl;
        out_complex_matrix(a.s,std::cout);
	std::cout << "va:" << std::endl;
        out_complex_matrix(a.v,std::cout);
	std::cout << "ta:" << std::endl;
        out_complex_matrix(a.t,std::cout);
        
        std::cout << std::endl;

	std::cout << "hb:" << std::endl;
        out_complex_matrix(b.h,std::cout);
	std::cout << "sb:" << std::endl;
        out_complex_matrix(b.s,std::cout);
	std::cout << "vb:" << std::endl;
        out_complex_matrix(b.v,std::cout);
	std::cout << "tb:" << std::endl;
        out_complex_matrix(b.t,std::cout);
        
//        m.random(9,7,7);

        std::cout << std::endl;

	std::cout << "hm:" << std::endl;
        out_complex_matrix(m.h,std::cout);
	std::cout << "sm:" << std::endl;
        out_complex_matrix(m.s,std::cout);
	std::cout << "vam:" << std::endl;
        out_complex_matrix(m.va,std::cout);
	std::cout << "tam:" << std::endl;
        out_complex_matrix(m.ta,std::cout);
	std::cout << "vbm:" << std::endl;
        out_complex_matrix(m.vb,std::cout);
	std::cout << "tbm:" << std::endl;
        out_complex_matrix(m.tb,std::cout);
#endif        


}


